% -----------------------------------------------------------------
% PARAMETRIC METHODS - AUTOREGRESSION
% Parametric_Methods.m
% MAR = Autoregression with moving window (fixed window size)
% AR = Autoregression with windows determined by the label
% Make sure the last row contains the labels
%------------------------------------------------------------------

  disp('Performing Autoregression');
    
  numSG=size(feature_matrix,1)-1; %number of features
 
  AROrder=10;
  AROffset = 1;    
  metricCount  = 0;
  metricOffset = 1; %AR
  %WindowSize = 20; %MAR
  %metricOffset = WindowSize; %MAR
 
   
  while(metricOffset<size(feature_matrix,2) )
    metricCount = metricCount + 1;
    
    %Normalize
    for i=1:numSG
    feature_matrix(i,:) = feature_matrix(i,:)-mean(feature_matrix(i,:));
    end;
    
    %determine window size (only in AR mode)      
    j=0;   
    k=metricOffset;
    WindowSize=1;        
    while (j == 0) && (k<size(feature_matrix,2))   
        k = k+1;
        if feature_matrix(numSG+1,k)==feature_matrix(numSG+1,k-1)
            WindowSize = WindowSize+1;            
        else
            j = 1;%end of window
        end;        
    end;
    
    % Perform Autoregression - BURG, COV, MCOV
    for i = 1:numSG        
      AROffset = (i*(AROrder+1))-AROrder;
 %     xx = feature_matrix(i,metricOffset-WindowSize+1:metricOffset); %image MAR
      xx = feature_matrix(i,metricOffset:metricOffset+WindowSize-1); %image AR
     
      %Determine AR coefficient
      [state.FV_burg(metricCount, AROffset:AROffset+AROrder), state.NV_burg(metricCount,i)] = arburg(xx, AROrder);
      [state.FV_cov(metricCount, AROffset:AROffset+AROrder), state.NV_cov(metricCount,i)] = arcov(xx, AROrder);
      [state.FV_mcov(metricCount, AROffset:AROffset+AROrder), state.NV_mcov(metricCount,i)] = armcov(xx, AROrder);
      
      %Determine AIC
      state.AIC_burg(metricCount, i) = WindowSize*log(state.NV_burg(metricCount,i))+2*AROrder;
      state.AIC_cov(metricCount, i) = WindowSize*log(state.NV_cov(metricCount,i))+2*AROrder;
      state.AIC_mcov(metricCount, i) = WindowSize*log(state.NV_mcov(metricCount,i))+2*AROrder;
        
    end; 
  
    % Determine entropy
    state.ENT_burg(metricCount) = entropy(state.FV_burg(metricCount,:));
    state.ENT_cov(metricCount) = entropy(state.FV_cov(metricCount,:));
    state.ENT_mcov(metricCount) = entropy(state.FV_mcov(metricCount,:));

    %add label
    state.FV_burg(metricCount,numSG*(AROrder+1)+1) = feature_matrix(numSG+1,metricOffset);
    state.FV_cov(metricCount,numSG*(AROrder+1)+1) = feature_matrix(numSG+1,metricOffset);
    state.FV_mcov(metricCount,numSG*(AROrder+1)+1) = feature_matrix(numSG+1,metricOffset); 
    
    metricOffset = metricOffset + WindowSize; %AR
   % metricOffset = metricOffset + 1; %MAR
        
  end % while
  
  %Remove first coefficient, since it is always 1
  i = numSG*(AROrder+1)-(AROrder);
  while i>0
       state.FV_burg(:,i)=[];
       state.FV_cov(:,i)=[];
       state.FV_mcov(:,i)=[];  
      i=i-(AROrder+1);
  end; 
  
  %Calculate mean AIC and entropy
   AIC_burg = mean2(state.AIC_burg);
   AIC_cov = mean2(state.AIC_cov);
   AIC_mcov = mean2(state.AIC_mcov);  
   ENT_burg = mean(state.ENT_burg);
   ENT_cov =mean(state.ENT_cov);
   ENT_mcov =mean(state.ENT_mcov);
    
  % Error Checking AR  
%   assert(metricOffset-1 == size(feature_matrix,2), 'ERROR: Processing of Feature Matrix Initial and Final NaN points fail');
%   assert(size(state.FV_burg,1) == (metricCount),'ERROR: Parametric Feature Matrix Dimensions invalid(row)');
%   assert(size(state.FV_burg,2) == (numSG*AROrder)+1,'ERROR: Parametric Feature Matrix Dimensions invalid(column)');
  
  clear  i j l xx k metricOffset AROffset metricCount AROrder numSG;